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SUMMARY 


The  emerging  field  of  nanofluidics  has  significance  to  the  Air  Force  with  application  to  nanopro¬ 
pulsion,  nano-rheology,  nano-lubrication,  control  of  nanosurface  properties  and  lab-on-a-chip  tech¬ 
nologies.  Biological  and  manufacturing  processes  that  involve  fluidic  transport  through  nanoscale 
membranes  and  channels  are  also  of  significance  to  the  Air  Force.  Design  and  operation  of  fluidic 
nanodevices  as  well  as  prediction  and  optimization  of  nanofluidic  processes  require  proper  mathemat¬ 
ical  and  computational  discretization.  Modeling  and  simulation  at  the  nanoscale  (mesoscopic)  range 
requires  new  mathematical  and  computational  methods  because  the  traditional  atomistic  (microscop¬ 
ic)  and  hydrodynamic  (macroscopic)  descriptions  are  not  valid  for  either  computational  or  theoretical 
reasons.  The  goals  of  research  performed  under  FA9550-06- 1-0236  were  to: 

•  Develop  the  theoretical  and  numerical  modeling  of  fluid  flow'  at  the  nanoscale  (mesoscopic) 
range. 

•  Address  new  multi  length-scale,  multi  time-scale  phenomena,  and  solid  wall-fluid  interactions 
that  arise  at  the  nanoscale. 

•  Investigate  fluidic  systems  that  incorporate  geometries,  fluids,  and  surface  materials  that  are 
representative  of  nanodevices  or  nanoprocesses  of  interest  to  the  Air  Force. 

The  Objectives  of  the  research  performed  under  FA9550-06-I-0236  were  to: 

•  Develop  advanced  particle  computational  methods  for  micro-  and  nano-scale  flows. 

•  Develop  numerical  boundary  conditions. 

•  Validate  the  numerical  methods  and  evaluate  their  numerical  errors. 

•  Investigate  flows  in  nanochannels  with  length  scales  from  10  nm  to  1000  nm. 

We  followed  two  distinct  but  inter-related  Technical  Approaches  under  FA9550-06-1  -0236: 

1)  Further  developed  the  unstructured  DSMC  (U3DSMC)  method  and  code,  originally  devel¬ 
oped  by  the  PI  under  a  previous  AFOSR  Grant  F49620-03-1-0219.  The  new'  computational 
developments  during  this  performance  period  allow: 

a)  high-fidelity,  gaseous  nanoscale  simulations; 

b)  evaluation  of  error  and  uncertainty  in  DSMC  nanoscale  simulations;  and, 

c)  comparisons  with  results  from  the  multiscale  Dissipative  Particle  Dynamics  (DPD)  simu¬ 
lations. 

2)  Developed  a  hierarchical,  multiscale,  top-down  approach  from  the  Smoothed  Particle  Hydro¬ 
dynamics  (SPH)  to  Smoothed  Dissipative  Particle  Dynamics  (SDPD)  that  allows: 

a)  computation  of  gaseous,  liquid,  and  complex  nanoflows; 

b)  a  two-level  (SDPD/SPH)  discretization  of  a  mesoscale  flow  system  by  changing  the  cha¬ 
racter  of  the  particle-particle  interactions  in  order  to  move  down  the  temporal-spatial 
scales. 
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The  Mathematical  and  Computational  Accomplishments  of  the  U3DSMC  approach  are: 

1 .  Explored  computational  issues  that  arise  in  the  DSMC  computation  of  nanoscale  flows. 

2.  Developed  the  mathematical  formulation  and  implemented  a  method  for  inflow  and  outflow 
boundary  conditions  in  U3DSMC  based  on  theory  of  characteristics. 

3.  Evaluated  the  numerical  uncertainty  and  statistical  fluctuations  in  macroscopic  flow  proper¬ 
ties  obtained  from  U3DSMC  computations  at  the  nanoscale.  The  statistical  error  in  number 
density,  mean  velocity  and  translational  temperature  obtained  from  U3DSMC  simulations 
was  compared  w  ith  theoretical  estimates.  The  investigation  considered  the  effects  of  the  num¬ 
ber  of  computational  particles  per  Delaunay  cell  and  Mach  number  on  the  fractional  error  for 
uniform  and  the  pressure  driven  nanoflows. 

4.  Performed  validation  and  grid  sensitivity  analysis  of  the  U3DSMC  method  at  the  nanoscale: 

a.  Compared  U3DSMC  simulations  of  supersonic  flows  into  microchannels  with  2D 
DSMC  simulations. 

b.  Compared  heat  transfer  rates  obtained  with  the  U3DSMC  method  at  the  upper  wall  of 
a  nanochannel  with  theoretical  estimates  applicable  to  the  free-molecular  regime. 

c.  Compared  mass  flux  at  the  outlet  of  a  nanochannel  with  the  results  from  the  semi- 
analytical  theory  developed  for  free-molecular  tube  flows. 

d.  Performed  a  grid  sensitivity  analysis  for  the  U3DSMC  method  with  simulations  in¬ 
volving  flows  into  nanochannels. 

5.  Applied  U3DSMC  to  explore  physical  phenomena  in  supersonic  sub-atmospheric  and  atmos¬ 
pheric  gas  flows  into  nanochannels  and  established  the  applicability  of  U3DSMC  at  the  na¬ 
noscale. 

The  Mathematical  and  Computational  Accomplishments  of  the  SDPD  approach  are: 

1.  Initiated  the  development  and  computational  implementation  of  a  hierarchical,  multiscale 
three-dimensional  computational  method  that  integrates  the  smooth  dissipative  particle  dy¬ 
namics  and  smooth  particle  hydrodynamic  techniques. 

2.  Developed  a  new  method  for  SDPD  solid-wall  boundary  conditions. 

3.  Performed  preliminary  validation  of  the  SDPD  method. 
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1  UNSTRUCTURED  3D  DIRECT  SIMULATION  MONTE  CARLO  (U3DSMC)  METHOD 


FOR  NANOSCALE  COMPUTATIONS 

We  review  below  the  elements  of  the  U3DSMC  method  developed  under  FA9550-06-1 -0236.  We 
summarize  also  the  numerical  challenges  associated  with  the  DSMC  simulations  of  nanoscale  physi¬ 
cal  domains.  Several  basic  algorithmic  modules  of  U3DSCM,  such  as  loading,  injection,  particle¬ 
tracing,  and  sampling  of  macroscopic  variables,  were  developed  by  the  PI  under  the  previous  AFOSR 
Grant  F49620-03-1-0219  and  are  also  shared  by  the  unstructured  particle  in  cell  a  methodology 
(U3dPIC)  developed  by  Gatsonis  and  Spirkin  (2009)’  for  the  simulation  of  fully  ionized  flows.  The 
new  developments  accomplished  during  3/20006-5/2009  encompass  the  mathematical  and  computa¬ 
tional  aspects  of  the  U3DSMC  particle  approach  for  high-fidelity  gaseous  nanoscale  simulations. 
These  developments  appear  in  Chamberlin  (2007)%  Al-Khouz  (2009)^  and  Gatsonis  and  Al-Khouz 
(2008)4.  The  U3DSMC  code  has  also  been  applied  to  the  simulation  of  gaseous  expansion  from  mi¬ 
cronozzles  (Chamberlin  and  Gatsonis,  2007)  and  the  simulation  of  microscale  jets  (Chamberlin  and 
Gatsonis,  2008)6. 

1 . 1  Computational  Issues  of  DSMC  for  Nanoscale  Flows 

The  increasing  proliferation  of  MEMS  devices  as  well  as  the  design  of  NEMS  devices  and  sen¬ 
sors  that  process  gases,  requires  modeling  and  analysis  of  flows  in  nanoscale  (or  submicron)  channels. 
The  physical  and  the  DSMC  computational  aspects  of  nanoflows  have  not  been  addressed  so  far  and 
have  some  unique  features.  The  computational  of  domains  of  interested  are  characterized  by  a  scale 
L  —  100  —  1000  nm  presents  two  challenges  to  the  implementation  of  the  DSMC  method. 

The  first  challenge  relates  to  the  small  number  of  real  particles  present  in  a  nanoscale  domain. 

Figure  1  depicts  the  real  number  of  particles  in  a  domain  with  volume  L 3  #  and  number  density  in  the 
range  of  n0  =  2.688  x  102j  m'1  (at  1  atm)  down  to  0.01  nQ.  The  number  of  particles  for  L  —  100  nm 
decrease  from  N  zz  104  at  nQ  to  N  «  102  at  0.01  nQ.  For  L  —  1000  nm  the  number  of  particles  be¬ 
comes  N  £3  1C)7  at  nQ  and  reduces  to  N  %  105  at  0.01nQ.  It  is  clear,  that  a  DSMC  computation  of  a 

nanoscale  domain  L3  may  require  particle  weights  that  are  much  smaller  when  compared  with  those 
used  traditionally  in  micro-  or  macroscale  applications,  resulting  in  even  fewer  computational  par¬ 
ticles  in  the  domain.  The  statistical  fluctuations  in  sampled  properties  inherent  to  microflows  is  ex¬ 
pected  to  be  more  important  to  nanoflows.2’  The  accuracy  of  the  computational  parameters  used  in 
the  nanoscale  simulations  (e.g,  particle  weights,  cell-sizing,  and  time-steps)  was  verified  with  a  grid 
sensitivity  study  and  comparisons  of  the  U3DSMC  results  with  theoretical  estimates. 
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Figure  1  The  Knudsen  number  as  function  of  nanoscale  length  L  and  the  number  of  real  particles  in  a  vo¬ 
lume  l) .  The  standard  number  density  /i0=2.69x102,;  m  \ 


The  second  challenge  relates  to  large  Knudsen  numbers  associated  with  nanoflows.  Figure  1 
shows  that  for  L  a  100  —  1000  nm  a  flow  is  expected  to  be  in  the  transitional  to  free-molecular  re¬ 
gime  for  all  densities  considered.  Even  under  atmospheric  pressure  conditions  (density  of  nQ)  the 
Kn  >  0.05.  While  transitional  to  near-free  molecular  conditions  are  often  encountered  in  rarefied 
microscale  or  macroscale  flows,  it  is  also  the  case  that  their  computational  cells  contain  a  large  num¬ 
ber  of  particles.  At  the  nanoscale,  the  entire  domain  of  interest  characterized  by  L  ,  can  be  of  the  or¬ 
der  or  smaller  than  A  .  A  DSMC  simulation  of  a  nanodomain,  therefore,  requires  discretization  with 
computational  cells  that  can  be  of  the  same  order  as  the  entire  domain  with  a  small  number  of  real 
particles  residing  into  it.  Furthermore,  the  number  of  sampling  cells  used  for  the  evaluation  of  ma¬ 
croscopic  fluid  variables,  such  as  density,  mean  velocity,  temperature,  may  not  be  sufficient  large  to 
provide  meaningful  spatial  resolution  in  a  nanodomain.  As  a  result,  flow  analysis  requires  the  use  of 
kinetic  properties,  such  as  the  phase-space  distributions,  rather  than  macroscopic  properties. 

1.2  Unstructured  3D  DSMC  Method 

The  U3DSMC  method  is  implemented  on  an  unstructured,  tetrahedral  Delaunay  mesh  (Chamberlin, 
2007).  The  grid  generator  discretizes  the  domain  fi  by  G0  Delaunay  tetrahedra  that  have  edge 
lengths  scaling  with  a  fraction  of  the  local  mean-free  path,  A .  Each  node  d  is  associated  with  the 
Voronoi  dual,  r ,  and  the  Delaunay  supercell  formed  by  all  the  Delaunay  cells  that  share  the 
node  d  as  illustrated  in  Figure  2  for  a  2d  geometry. 
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Figure  2.  Grid  structure  used  in  the  U3DSMC  depicted  in  2D  for  clarity.  The  Delaunay  cell  11^,  the  Vo- 
ronoi  dual  ,  and  a  Delaunay  supercell  and  a  node  d . 

1.2.1  Particle  L  oading  and  Injection 

A  number  of  computational  particles,  each  representing  F  real  particles,  is  loaded  in  the  domain  and 

injected  from  boundaries.  In  U3DSMC,  particle  loading  is  carried  out  by  placing  species  s  particles 
in  each  Delaunay  cell  with  randomly  chosen  positions  and  velocities  following  the  quasi-equilibrium, 
drifting  Maxwellian  distribution 

\3/2 


2t rfcT 


exp 


w.(v  -V) 


2  kT 


(1) 


where  number  density  ris(r.t),  translational  temperature  T  (r A),  and  mean  species  velocity  V  (r, t) 
are  given.  The  number  of  computational  particles  loaded  at  each  Delaunay  cell  is 

ND=n(r,t)QD/Fw  ^ 

where  Fw  is  the  particle  weight.  Computational  particles  are  injected  inside  the  domain  from  its 

boundaries.  The  inward  flux  of  species  s  real  particles  across  a  Delaunay  boundary  surface  is  due  to 
a  drifting  Maxwellian  eq.  (1) 

nCm  ( 


N 


~LJ^=~  (exP  S/cos2  (#))  +  >/tt  S9  cos  (0) |l  +  erf  cos  Jjj 


(3) 


where  the  most  probable  thermal  speed  is 

C,ra  =  PkHTJm,  (4) 

the  mean  velocity  V  is  inclined  at  an  angle  0  to  the  unit  normal  vector  e .  The  speed  ratio  of  the 
mean  speed  to  the  most  probable  speed  is 

S.  =  V/C :  (5) 
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The  number  of  species  s  computational  particles  AN*  to  be  injected  to  the  domain  in  a  given  time 
step,  is 


(6) 


where  Ar  is  the  time  step  and  A  is  the  area  of  the  surface  element. 

1 .2.2  Integration  of  Particle  Motion  and  Search  for  Particle  Location 

The  computational  particles  are  tracked  as  they  interact  with  other  particles  and  with  domain 
boundaries.  For  a  particle  p  with  position  rp(t)  and  velocity  vp(£)  =  (v  j,t?  ,  t;  )  the  new  particle 
position  is  obtained  from  integration  of  the  equation  of  motion,  given  by 
r  (t  +  dt)  =  rp(t)  +  vp  (j)  At.  Tracing  the  motion  of  the  particle  on  the  unstructured  grid  is  carried 

out  using  the  successive-neighbor  search  algorithm  (Gatsonis  and  Spirkin,  2008) 1 .  In  U3DSMC  mole¬ 
cular  motion  over  the  time  step  (Ar)  and  intermolecular  collisions  occurring  at  a  mean  collision 
time,  r  ,  are  uncoupled,  by  choosing  Ar  <  r . 

1.2.3  Elastic  and  Inelastic  Collisions 

The  selection  of  collision  partners  is  performed  within  each  Delaunay  cell  and  follows  Bird’s  no¬ 
time-counter  (NTC)  scheme/  The  cross-section  for  elastic  collisions  is  based  on  the  variable  hard 
sphere  model  (VHS).  The  Larsen-Borgnakke  model  is  used  to  simulate  the  exchange  of  rotational 
energy  between  the  collision  pair  Solid  surfaces  are  modeled  as  being  diffuse,  non-diffuse  or  spe¬ 
cularly  reflecting.  (Chamberlin,  2007)2 

1 .2.4  Inflow  and  Outflow  Boundary  Conditions 

Inlet  and  outlet  boundary  conditions  implemented  in  U3DSMC  can  accommodate  supersonic  and 
subsonic  flows.  For  supersonic  flows  particles  are  injected  using  eq.  (6)  based  on  specified  n  (r,£), 

T  (r,  t),  and  V(r ,t).  The  subsonic  boundary  method  implemented  in  U3DSMC  follows  an  ap¬ 
proach  that  is  based  on  the  method  of  characteristics  used  in  compressible  flow  computations.  Similar 
methods  have  been  used  in  2D  DSMC  subsonic  simulations  (Nance  et  al.,  1997; 10  Liou  and  Fang, 
2000;11  Wang  and  Li,  2004 12).  Our  method  was  generalized  to  3d  unstructured  grids  and  was  sup¬ 
plemented  with  a  neighboring-cell  sampling  approach  and  a  smoothing  technique  to  speed  conver¬ 
gence.  We  followed  Flirsch  (1995)  and  Yee  at  al.  (1982) 14  and  derived  the  boundary  methodology 

implemented  in  U3DSMC.  At  a  subsonic  inflow  boundary,  the  characteristic  variables  wvtv2 

represent  the  physical  boundary  conditions,  while  at  a  subsonic  outflow,  w  is  the  physical  boundary 
condition.  For  a  subsonic  inlet  the  two  physical  boundary  conditions  specified  at  the  boundary  are 
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and  Pj .  Zero-order  space  and  zero-order  time  extrapolation  of  the  characteristic  variable  w3  be¬ 
tween  the  boundary  point  and  the  interior  point  provides  the  numerical  boundary  condition  and  allows 


the  evaluation  of  the  third  primitive  variable,  Vx .  For  a  subsonic  outlet  designated  by  subscript  M  , 


w3  is  the  physical  boundary  condition  and  the  pressure  is  fixed  at  pM .  Zero-order  space  and  zero- 
order  time  extrapolation  of  the  characteristic  variables  and  use  of  the  numerical  boundary  conditions, 
provide  the  two  primitive  variables  pM  and  V  .  Once  n,  T,  V  are  obtained  on  an  inlet  or  outlet  sur¬ 
face,  the  particles  are  injected  into  the  domain  following  the  injection  procedure  of  eq.  (6). 

1.2.5  Sampling  and  Macroscopic  Flow  Properties 

In  U3DSMC  macroscopic,  single-fluid  and  multi-fluid  variables,  such  as  density,  mean  velocity, 
etc,  are  evaluated  on  nodes  by  averaging  the  particle  properties  within  the  Delaunay  supercell  or 

the  Voronoi  cell  Td  shown  in  Figure  2.  For  A faD{t)  computational  particles  of  species  s  residing  in 
the  sampling  Delaunay  cell  with  volume  the  number  density  is  defined  as 


(7) 


the  mean  species  velocity  of  the  cell  V  (Z),  £)  =  V^,  is 


(8) 


and  the  mass-average  velocity  of  the  cell  V(A  t)  —  jlZ,  Vy)  is 


(9) 


The  thermal  velocity  Csp  is  defined  with  respect  to  the  mass-average  velocity  V  as 


C  =v  -V 


(10) 


The  translational  temperature  is 


P=i 


(11) 


and  the  scalar  pressure  is 

p9{D9t)  =  n(D4)kT{D,t). 


(12) 


We  also  calculate  the  overall  single-fluid  cell  properties  as 


T(D,t)  =  '£/nT;(D,t)/'£n(D,t) 


(13) 
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p(D,t)  =  ^Pl 


(14) 


In  addition  to  such  cell-averaged  properties  the  instantaneous  value  of  any  macroscopic  variable 
X(d,t)  =  n(dit),  V(d,i),Trf(d,i)etc  associated  with  a  node  d  shown  in  Figure  3,  is  obtained  with  a 

volume  averaged  procedure  in  the  Delaunay  supercell  Q ,  following, 


CCH3  j 

x(d,t)  =  J2x(Dj)nD  nd 


(15) 


For  steady  U3DSMC  simulations  we  obtain  after  reaching  steady-state  m  =  1,  M  independent  cell- 
based  or  nodal -based  samples.  The  sample-averaged  mean  macroscopic  property  for  a  node  (or  cell) 
is  denoted  by  X^D) 

M 


rUoy v,«dvt^d)  etc and is siven by 


=  2. 

m=l  / 


(16) 


Several  flux  properties,  such  as  pressure,  shear  stress  and  heat  flux  (W/m2)  to  a  solid  surface,  are 
determined  from  the  perpendicular  or  parallel  momentum  and  total  energy  exchange  of  the  impinging 
particles.  For  example  for  heat  flux, 


A  At 


(17) 


In  the  above,  E  is  the  reflected  translational  and  internal  energy  of  impinging  particles  on  surface 
area  A  ,  and  At  is  the  duration  of  impingement  sampling. 


1.3  Validation  of  U3DSMC  for  Micro-scale  Flows 

The  U3DSMC  code  was  validated  for  micron-scale  computations  by  comparisons  with  previous  2d 
DSMC  results  of  Tiou  and  Fang  (2001) 15  and  Le  et  al.  (2007)16  These  simulations  involved  super¬ 
sonic  nitrogen  flow  incoming  into  a  microchannel  with  height  of  1.2  pm,  length  6  pm  and  represent 
the  smallest  scales  available  in  the  literature  for  comparison  The  domain  used  in  the  U3DSMC  simu¬ 
lations  is  depicted  in  Figure  2.  The  incoming  nitrogen  N2  flow  has  nM  =  1.75  x  1025  m‘\ 

T  =  300K  ,  V  =  1465.7  m/s  corresponding  to  P  =  72450  Pa,  M  —  4.15,  and  A  =  74.4  nm, 

Kno c  =0.062.  Following  [15,  16]  the  first  0.1  L  of  the  length  of  the  computational  domain  is  as¬ 
sumed  to  be  specularly  reflecting  in  order  to  generate  uniform  free-stream  condition  inside  the  micro- 
channel.  The  supersonic  incoming  conditions  are  applied  at  the  inlet  of  the  computational  domain.  A 
vacuum  outlet  boundary  condition  is  specified  at  the  exit  of  the  microchannel.  Particles  that  leave  the 
exit  are  removed  from  the  domain  and  no  particles  are  allowed  to  return  into  the  domain.  The  upper 
and  lower  walls  of  the  microchannel  are  isothermal  at  a  constant  temperature  of  323  K  and  are  as¬ 
sumed  to  be  fully  diffuse  with  perfect  thermal  accommodation.  The  side  walls  are  simulated  as  specu- 
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lar  so  that  they  act  as  symmetry  planes.  In  order  to  compare  with  the  2  D  results,  the  depth  D  of  the 
microchannel  in  our  3D  simulations  was  varied  and  it  was  found  that  results  are  not  sensitive  for 

D  >  A  .  A  D  =  74.4  nm  was  used  in  the  simulations.  The  domain  is  discretized  with  23032  Delau- 

-  oc 

nay  cells  with  405  particles/cell  upon  loading  and  At  =  1x10  13  sec.  Macroscopic  flow  properties 
are  sampled  once  steady  state  is  reached.  Twenty  samples  are  used  in  obtaining  the  average  after 
reaching  the  steady  state  time. 

< 

W 

A 
A 

H  Exlt 

Region 

T 

◄  x  .  _ ,  ► 

0.1  L  0.9  L 

Figure  3.  Geometry  used  in  the  U3DSMC  simulation  of  a  microchannel  flow  for  comparison  with  2D 
DSMC  simulations  of  Liou  and  Fang  15  and  Le  et  al.  lft 

Figure  4  shows  the  centerline  y  /  H  =  0.5  translational  temperature  profile  given  in  Eq.  (11)  obtained 
from  the  U3DSMC  simulation.  The  difference  between  Liou  and  Fang,  Le  et  al.  and  our  3d  results  is 
less  than  7%.  Figure  4b  shows  an  excellent  agreement  between  the  velocity  profiles  obtained  for 
x  /  L  =  0.5  .  Figure  4c  shows  a  good  agreement  with  the  heat  flux  at  the  wall.  Differences  in  heat 

flux  rates  between  the  3d  and  2d  simulations  can  be  attributed  to  the  sidewalls  present  in  our  3d  simu¬ 
lations  as  well  as  the  sampling  techniques  used  to  obtain  this  surface  fluxal  property. 

1.4  Statistical  Error  and  Fluctuations  in  the  U3DSMC  Method 

In  a  dilute  gas  there  are  inherent  statistical  fluctuations  in  macroscopic  parameters  such  as  number 
density,  mean  velocity,  temperature,  that  are  provided  by  theoretical  estimates.  In  the  case  of  a  par¬ 
ticle  simulation  such  as  the  U3DSMC,  these  statistical  fluctuations  result  to  statistical  errors  due  to  the 
small  number  of  particles  used  to  represent  the  real  number  of  particles  in  the  domain,  and  the  sam¬ 
pling  associated  with  estimation  of  macroscopic  variables.  Such  fluctuations  and  errors  when  com¬ 
pared  to  theoretical  estimates,  are  functions  of  the  particle  weight  ( Fw ),  the  volume  of  the  Delaunay 

cells  the  time-step  (  At  )  used,  and  the  number  of  samples  ( M )  used  to  obtain  macroscopic 

properties.  We  followed  Hadjiconstantinou  et  al.  (2000) 1  and  generalized  their  derivation  for  3D  un¬ 
structured  grids  used  in  U3DSMC  considering  a  gas  of  species  s  particles. 
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Figure  4.  Temperature  along  centerline,  mean  velocity  at  a/L=0.5,  and  heat  flux  to  the  wall,  in  a  micro- 
channel  flow.  Comparison  between  U3DSMC  with  2D  DSMC  results  of  Liou  and  Fang  [15|  and  Le  et 
a/.  [161 
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We  considered  the  fluctuations  associated  with  three  important  macroscopic  fields,  the  number 
density  ( n ),  mean  velocity  ( V  )  and  translational  temperature  ( T  )  obtained  over  m  =  1,  M  samples 
as  shown  in  eq.  (16)  The  sample  standard  deviation  for  any  cell-based  variable  at  the  m- th  sample 
X™  =  n™ ,  V"‘ ,  T™  is  defined  as 


i 


M  -  1 


(18) 


and  the  standard  deviation  (uncertainty)  of  any  sample-averaged  mean  variable  XD  =  nD,\D,TD  ob¬ 
tained  from  eq.  (16)  is 

<r(XD)  =  a(X;)/X\  (19) 


10 


The  fractional  error  of  a  sample-averaged  mean  variable,  is  defined  as  the  standard  deviation  of  the 
mean  over  the  mean 


E 

x°  xD 

The  theoretical  fractional  error  in  the  number  density  in  a  Delaunay  cell  is 
E  =  1  1 


(20) 


(21) 


"  no  nd  ac 

In  the  above  equation  the  acoustic  number  is  Ac  =  aj a .  ,  where  a  is  the  speed  of  sound,  and 


m  is  the  speed  of  sound  at  a  reference  temperature.  For  dilute  gases  Ac-\.  Similarly, 
the  fractional  error  in  the  mean  velocity  component  V  is 


E„  = 


(22) 


^/maT  AcMyfi 

where  1  =  cp  /  cv  is  the  ratio  of  specific  heats  and  the  Mach  number  is  M  =  V^/a.  The  fractional 
error  in  temperature  is 


"fe) 


Er  = 

T. 


(23) 


tJ m  ,/maTVc. 

where,  C  is  the  heat  capacity  at  constant  volume.  The  equations  (21 ),  (22),  and  (23)  provide  also  the 
statistical  error  for  sample-averaged  nodal  variables  Xd  =  ftp  ,7^  based  on  M  instantaneous  nod¬ 
al  samples,  X™  =  n™  .  The  same  equations  provide  also  the  fractional  error  in  the  real  gas 


if  used  with  the  real  number  of  particles  in  the  Delaunay  or  the  super-Delaunay  cell,  =  ND^FW . 


1 .4. 1  Supersonic  and  Subsonic  Nanoscale  Flows 

In  order  to  evaluate  the  statistical  fractional  errors  in  U3DSMC  nanoscale  simulations,  we  per¬ 
formed  simulations  for  a  flow  of  N2  with  n  =2.69xl025  rn"3,  T  =  273  K,  V  =  3369.17  m/s, 

=  10  with  48.1  run  ,  in  a  domain  that  requires  the  discretization  of  a  characteristic  length 

scale  L  =  50 nm.  This  set  of  conditions  corresponds  to  one  of  the  most  restrictive  for  U3DSMC 
computation  cases  due  to  the  smallest  mean-free  path  involved.  These  U3DSMC  simulations  were 
realized  in  a  rectangular  domain  of  (H  xW  x  D)  region  with  II  =  W  =  0.05  pm  an  dD  =  10  pm . 

The  simulations  were  run  with  At  =  1  x  10  12  s  and  steady  state  was  achieved  after  104  time  steps 
by  monitoring  the  mass  flow  rate  at  the  exit.  The  domain  was  discretized  with  GD  =  2446  Delaunay 
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cells  with  an  average  volume  (p>D}  =  9.9x10  24  m3,  average  edge-length  =  0.0438 //in and  an 

average  number  of  real  particles  per  cell,  ND  —  266  .  Two  weights  were  used  in  the  simulations, 

Fw  —  26  with  N D  =  10  computational  particle  per  cell,  and  Fw  =  9  with  N  D  =  30.  The  results 
from  the  U3DSMC  simulations  provided  M=  460  Delaunay  samples  that  were  used  to  obtain  sam¬ 
ple-averaged  nD,VD, T0  following  eq.  (16).  The  standard-deviation  a x  is  evaluated  from  eq.  (19), 

and  subsequently  eq.  (20)  provides  the  numerical  error,  designated  as  Ex  (U3DSMC) .  We  also  eva¬ 
luated  the  fractional  error  using  eq.  (21),  (22),  (23)  with  U3DSMC  parameters,  and  designate  as 
£X(U3DSMC  Analytical) .  Finally,  we  evaluated  the  error  /^(Real  Analytical)  following  equa¬ 
tions  (21),  (22),  (23)  using  the  real  number  of  particles  in  the  sampling  cell  . 

Figure  5  shows  the  fractional  error  as  a  function  of  the  computational  particles  N D  and  thus  par¬ 
ticle  weight.  The  results  show  that  U3DSMC  simulations  using  particle  weights  as  small  as  9,  have 
numerical  fractional  errors  that  compare  very  well  with  the  analytical  U3DSMC  errors.  It  is  therefore, 
desirable  to  use  the  theoretical  formulas  eq.  (21),  (22),  (23)  in  estimating  the  fractional  errors  in  a 
U3DSCM  simulation  instead  of  calculating  the  computationally  expensive  numerical  fractional  er¬ 
rors.  The  results  from  Figure  5  also  show'  that  the  theoretical  fractional  errors  in  the  real  gas  are  with¬ 
in  the  same  order  of  magnitude  as  the  numerical  errors.  These  results  show  that  the  no-time-counter 
(NTC)  collision  methodology  implemented  in  U3DSCM  is  robust  even  with  such  small  particle 
weights.  If  the  particle  weight  becomes  1 ,  then  certain  algorithmic  steps  in  NTC  could  require  mod¬ 
ifications  that  are  outside  the  scope  of  this  work  Figure  6  illustrates  the  effect  of  Mach  number  on 
the  fractional  error  Ev  using  the  same  number  of  samples  as  before.  The  fractional  error  for  the  sub¬ 
sonic  case  M  =  0.1  is  two  order  of  magnitudes  larger  than  the  supersonic  M  =  10  a  direct  result  of 
eq.  (22).  Our  results  are  consistent  with  similar  findings  in  structured  DSMC  computations  found  in 
Hadjiconstantinou  et  al.  (2003)17. 

1.4.2  Subsonic  Nanoflow 

A  pressure  driven  subsonic  flow  of  N2  was  also  simulated  using  the  U3DSMC.  The  geometry 
that  is  considered  is  a  nanochannel  of  500  nm  height  and  100  nm  width  4  //m  length  We  assumed  a 

mean  flow  with  =  2.69  xlO25  m  3,  Tx  =  273  K  corresponding  to  atmospheric  conditions  with 

a  =  101325  Pa  (1  atm)  and  48.1  nm.  WithPn,  and  Tm  specified,  Vm  was  obtained  from 

the  inflow  boundary  condition  method  described  in  Sec.  1.2.4,  and  used  for  particle  injection  at  the  in¬ 
let  boundary  during  the  iteration 
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Figure  5.  Statistical  fractional  errors  in  sample-averaged  macroscopic  density,  velocity  and  temperature 
In  a  nanoseale  flow.  The  U3DSMC  errors  are  a  function  of  the  number  of  computational  particles  in  a 
Delaunay  cell.  The  real  analytical  error  is  based  on  the  real  number  of  particles  in  a  Delaunay  cell. 
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Figure  6.  Statistical  error  in  mean  velocity  for  the  uniform  flow  as  a  function  of  the  Maeh  number. 
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Figure  7  shows  the  fractional  error  in  the  number  density,  mean  axial  velocity,  and  temperature 
as  a  function  of  the  transverse  direction  ( y  ).  The  U3DSMC  errors  are  almost  identical  with  the  errors 
obtained  from  the  theoretical  expressions  using  the  local  values  obtained  from  the  U3DSMC  and  the 
real  number  of  particles.  The  error  in  density  and  temperature  shows  no  dependence  on  the  transverse 
position.  The  error  in  mean  velocity  depends  on  the  position.  The  minimum  error  occurs  at  the  center 
of  the  channel  where  the  maximum  velocity  is  achieved 
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Figure  7.  Statistical  error  in  the  U3DSMC  simulation  of  a  pressure  driven  nanoflow. 
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1 .5  U3DSMC  Computation  of  Nanoscale  Flows 


Gas  flows  in  microscale  channels  have  been  studied  extensively  due  to  their  relevance  to  micro- 
machined  and  MEMS  devices  or  sensors.  In  the  case  of  continuum  and  slip  flow  regimes  characte¬ 
rized  by  the  inlet  Knudsen  number  based  on  the  height  H  of  the  microchannel  Kn^  =  \x  j H ,  analy¬ 
sis  and  simulations  are  usually  carried  out  with  Navier-Stokes  type  of  solvers.  In  the  case  of 
transitional  and  free-molecular  regimes  (  Kn >  0.01 ),  analysis  and  simulations  are  carried  out  with 

Boltzmann  solvers,  most  often  based  on  the  Direct  Simulation  Monte  Carlo  (DSMC)  method.  The 
DSMC  has  undergone  various  algorithmic  improvements  (Bird  et  al.,  2009)18  since  its  introduction 
and  has  become  the  primary  method  for  simulating  gaseous  flows  in  microdomains.(Karniadakis  et 
al.,  2005 ,9;  Liou  and  Fang,  200620.  While  there  is  a  large  number  of  investigations  covering  subsonic 
flows  in  microchannels,  there  is  a  smaller  number  covering  supersonic  microflows  and  very  few  in¬ 
vestigations  cover  the  sub-micron  or  nanoscale  regime.  The  research  performed  under  this  grant  pro¬ 
vided  a  First  insight  on  supersonic  flows  in  nanoscale  rectangular  channels.  The  simulations  per¬ 
formed  with  the  U3DSMC  code,  covered  the  transitional  to  near-free  molecular  regimes.  (Al-Khouz, 
2009)3. 

Rarefied  supersonic  flows  through  channels  and  tubes  has  been  studied  analytically2  and  compu¬ 
tationally  using  DSMC  due  to  its  importance  in  many  technical  applications.2"  23' 24  For  microchan¬ 
nels  the  2d  DSMC  simulations  so  far  have  addressed  the  slip  to  transitional  regimes  for  nitrogen  and 
helium  gases.  2  '  ^  2  '  16  MicroChannel  heights  considered  in  these  investigations  are  in  the  range 

0  1-2.4  microns  with  aspect  ratios  of  2.5-5  and  Knudsen  numbers  in  the  range  0.0046  to  0.744  cover¬ 
ing  the  slip  and  translational  regimes.  These  DSMC  investigations  so  far  considered  fully  diffuse 
walls.  In  Mavriplis  et  al.(1997)26  and  Oh  et  al.  ( 1 997)2  the  wall  temperature  was  specified  to  be 
equal  to  the  free  stream.  In  Liou  and  Fang  (2001  )15,  Le  and  Hassan,  (2006)2  and  Le  et  al.  (2007)16 
the  wall  temperature  was  set  higher  than  the  free  stream  for  parts  of  the  microchannels  to  investigate 
the  heat  transfer  phenomena.  These  investigations  have  shown  that  a  rarefied  supersonic  flow  incom¬ 
ing  into  a  macro-  or  microchannel  with  height  H  and  length  L  is  influenced  by  the  free-stream 

Knudsen  number  Kn ^  ,  the  free-stream  ratio  of  mean  speed  over  the  most  probable  speed  5  ,  the  as¬ 
pect  ratio  of  the  channel  All  =  L  /  //  ,  the  wall  temperature,  the  nature  of  wall  reflection,  and  the 
type  of  molecules  involved. 


1.5.1  Free  stream ,  Boundary  Conditions  and  Computational  Parameters 

The  supersonic  flow  into  a  rectangular  nanochannel  was  investigated  with  the  U3DSMC  code  for  the 
geometry  shown  in  Figure  8.  The  simulation  domain  featured  a  buffer  region  with  length  Lgy  shown 
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in  Figure  8,  followed  by  the  nanochannel  with  length,  height  and  width  L,  H.  W  respectively.  The 
flow  is  characterized  by  the  Knudsen  number 


Kn  =  —  = 
H 


1 


^d,xT.jT) 


W  i/2  // 


(24) 


with  d^f  —  4.17  x  10  10  m,  Tref  —  273  K  and  uj  =  0.74  181  The  local  speed-ratio  is  defined  in  eq.  (5) 
and  the  local  Mach  number  is  given  in  terms  of  the  specific  heat  ratio  7  =  1.4  ,  particle  mass 
m  =  4.68  x  10  26  kg,  as 


M  =  V 


7fcr 


m 


(25) 
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Figure  8.  The  computational  domain  used  in  the  U3DSMC  simulations  of  nanochannel  flows  (dimensions 

not  to  scale). 

Rectangular  nanochannels  with  H  —  W  =  100  —  1000  nm,  AR=1,  10,  100  are  considered  in  the  si¬ 
mulations  with  input  parameters  shown  in  Table  1.  The  incoming  nitrogen  N2  flow  has 

nx  =  2.69 x  1024rn  3 ,  Tx  =  273K,  corresponding  to  Px  =  O.latm  (10.132  kPa)  and  Ax  =  481 
nm.  Cases  1-9  have  V  =2013  m/s,  S  =5  and  M  =5.97.  Case  10  has  V  =  805.2  m/s,  S  = 

OC  X  X  x  x 

2,  M  =2.39.  Case  11  has  V  =  4026.8  m/s,  S  =10  and  M  =  11.95.  In  the  simulations  we  con- 

OC  (X  X  X 

sider  the  rotational  but  neglect  the  vibrational  degrees  of  freedom  of  the  nitrogen. 

In  order  to  allow  for  the  undisturbed  free  stream  conditions  to  be  realized  far  from  the  nanochan¬ 
nel  inlet  we  included  a  buffer  region  of  length  LB  ahead  of  the  nanochannel  inlet.  Molecules  were  in¬ 
jected  through  the  boundaries  of  the  buffer  region  based  on  eq.  (6)  with  specified  ,  7^ ,  .  The 
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sides  of  the  buffer  region  were  modeled  as  a  free  stream.  Particles  that  moved  upstream  and  reachEd 
the  surfaces  of  he  buffer  region  were  removed  from  the  computational  domain. 


Table  1.  Physical  and  computational  parameters  for  U3DSMC  simulations  of  nitrogen  subatmospheric 
nanochannel  flows. 


Case 

H 

(nm) 

LfH 

s 

oc 

M 

OC 

p> 

(kPa) 

Kn 

OC 

Lb 

(nm) 

Gv 

(/„)(nm) 

K) 

t»K)] 

1 

100 

1 

5 

5.97 

0 

4.81 

1000 

132 

85.8 

20 

11.21 

2 

100 

10 

5 

5.97 

0 

4.81 

1000 

252 

[8.74] 

85.2 

20 

[4.44] 

10.67 

3 

100 

100 

5 

5.97 

0 

4.81 

1000 

504 

[7.2] 

122 

20 

[3.78] 

29.36 

4 

500 

1 

5 

5.97 

0 

0.962 

1000 

1544 

[17.8] 

118 

20 

[10.96)] 

32.67 

4a 

500 

1 

5 

5.97 

120 

0.962 

100 

1544 

[12.1] 

118 

20 

[11.93] 

32.67 

4b 

500 

1 

5 

5.97 

200 

0.962 

100 

1544 

[12.1] 

118 

20 

[11.93] 

32.67 

5 

500 

10 

5 

5.97 

0 

0.962 

1000 

5836 

[12.1] 

129 

20 

[11.93] 

34.57 

5a 

500 

10 

5 

5.97 

40 

0.962 

1000 

5836 

[11.4] 

129 

20 

[11.03] 

34.57 

5b 

500 

10 

5 

5.97 

100 

0.962 

1000 

5836 

[11.4] 

129 

20 

[11.03] 

34.57 

6 

500 

100 

5 

5.97 

0 

0.962 

1000 

48576 

[11.4] 

135 

20 

[11.03] 

35.30 

7 

1000 

1 

5 

5.97 

0 
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The  outlet  of  the  nanochannel  was  placed  at  the  exit  boundary  of  the  computational  domain  as 
shown  in  Figure  8.  A  vacuum  boundary  condition  was  set  for  Cases  1-11,  where  the  back  pressure  is 

specified  as  Pb  =  0  at  the  exit  region.  Any  particle  reaching  the  exit  plane  was  removed  from  the  do¬ 
main  and  no  particles  were  allowed  to  enter  the  computational  domain  through  the  exit  plane.  For 
Cases  4a,  4b  and  5a,  5b  a  finite  back  pressure  Pb  was  specified  at  the  exit.  The  upper  and  lower  walls, 
as  well  as  the  side  walls  of  the  nanochannel,  were  modeled  as  fully  diffuse  and  the  temperature  of  the 


17 


wall  equals  that  of  the  free  stream  distribution.  The  simulation  domain  was  loaded  initially  with  a 
uniform  drifting  gas  following  eq.  (6)  based  on  the  free  stream  conditions  in  Table  1 . 


Figure  9.  Typical  grid  used  in  U3DSMC  simulations  (Case  4)  showing  the  surface  elements  of  the  Delau¬ 
nay  tethrahedra. 

Figure  9  displays  the  boundary  faces  of  the  Delaunay  cells  of  a  typical  mesh  used  in  the  simulations. 
The  interior  is  discretized  with  Delaunay  tethrahedra  cells  with  edge-lengths  that  are  smaller  than  A^  . 
The  mean  of  the  Delaunay  edge  lengths  (lD^  in  Table  1  show  that  they  are  about  one  fourth  of  the 

corresponding  Ax  —  481  nm.  The  total  number  of  Delaunay  cells  in  the  domain  G^,  the  average 

number  of  simulated  particles  in  each  cell  and  its  standard  deviation  cr(ND)  are  provided  in 

Table  1 .  These  parameters  show  that  the  Delaunay  cells  are  populated  at  initialization  with  at  least  10 
computational  particles. 

The  collision  time  under  free  stream  conditions  is  on  the  order  of  10  9  s  and  can  become  on  the 

order  of  10” 11  s  for  higher  densities  and  temperatures  encountered  in  the  nanochannels  considered  in 
these  simulations.  From  the  parameters  in  Table  1  the  most  restrictive  convective  time  scale 

{Id)  /  Kc  ls  on  order  °f  10  11  sec.  The  simulations  were  run  with  a  time  step  At  —  10“12  s  to  en¬ 
sure  decoupling  of  particle  motion  with  collisions,  and  ensure  that  computational  particles  do  not 
cross  a  cell  prior  to  a  collision.  The  steady  state  time  was  determined  from  the  mass  flow  rate  ob¬ 
tained  at  the  outlet  of  the  nanochannel,  as  shown  in  Figure  10.  After  reaching  the  steady  state  time  of 
0.2  psec,  100  nodal  samples  obtained  following  the  volume  average  procedure  of  eq.  (15),  corres¬ 
ponding  to  more  than  2,000  Delaunay  samples,  were  averaged  to  produce  the  steady-state  flow  prop¬ 
erties  following  eq.  (16).  The  10  11  seconds  (100  time  steps)  allowed  between  two  successive  nodal 
samples  combined  with  the  large  number  of  Delaunay  samples  used  in  obtaining  each  nodal  sample, 
provides  some  assurance  for  the  statistical  independence  of  the  samples. 
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Figure  10.  Mass  flow  rate  from  U3DSMC  simulations  at  the  outlet  of  a  rectangular  nanochannel  with 
—  0.1  atm  ,  aspect  ratio  AR,  and  height  H  (Case  1-0).  The  dotted  line  shows  the  steady-state  time 
used  to  obtained  samples  for  the  evaluation  of  macroscopic  variables. 


1.5.2  Grid  Sensitivity  A  naiysis 

The  cell  size  is  an  important  parameter  in  a  DSMC  simulation  and  for  microflows  it  has  been  ex¬ 
amined  extensively  for  2D  uniform  grids. 2K‘  1  For  the  unstructured  3D  simulations  in  this  work  the 
grid  was  generated  using  nearly  equally  sized  tetrahedra.  While  the  distribution  of  edge-lengths  is 
nearly  uniform,  the  mean  (lD^  is  the  proper  representative  cell  scale.  The  grid  sensitivity  analysis  stu¬ 
died  the  effect  of  increasing  the  cell  size  using  as  baseline  parameters  presented  in  Table  1  and  com¬ 
paring  macroscopic  flow  parameters  such  as  density,  flow  velocity  and  translational  temperature. 
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Table  2.  U3DSMC  simulation  parameters  used  in  the  grid  sensitivity  study. 


Case 

CT(/C)(nm) 

K) 

°(Nd) 

2 

85.2  [1/5.6] 

8.74 

11.21 

4.44 

GS-2 

136  [1/3.4] 

1 1.7 

35.87 

13.73 

4 

129  [1/3.4] 

11.4 

34.57 

11.03 

GS-4 

205  [1/2.3] 

17 

57.5 

11.2 

Results  for  the  two  most  restrictive  Cases  2  and  4,  i.e.  nanochannels  with  the  small  inlets  and  aspect 
ratios,  are  discussed  bellow.  Input  parameters  of  Case  2  and  Case  4  are  used  in  simulations  GS-2  and 

GS-4  with  nearly  doubling  the  edge-length  lD  while  keeping  the  particle  weight  constant  at  Fw  =  20  . 
The  resulting  (lD^  and  are  shown  in  Table  2. 


V(rrv*) 

7(K) 
n(' n'9) 

|GS) 

T  (K)  (G$) 
n  (m'*)  (GS) 


Case  2 


iaoo 

31000 


*1400 

1200 

1000 

MO 


- 

1500 

1400 

\j\  /  \ 

\\l  \ 

1200 

A  v\  ■ 

r\  i  i 

J  \  \ 

1Q00„ 

MO  ^ 

MO 

400 

200 

4e*25£ 


OA  10  1.9 

Case  4 


5**24 

7**24 

5**24 

WT 

5e*24jE 

c 

4**24 

3**24 

2**24 


*(H™) 

Figure  11.  Centerline  values  of  n,  V  ,T  .  Top:  Case  2  (lines)  and  Case  GS-2  (symbols)  Bottom:  Case  4 

(lines).  Case  GS-4  (symbols). 


Centerline  values  for  the  number  density,  translational  temperature,  and  mean  velocity  are  shown 
in  Figure  11.  Case  2  and  Case  GS-2  represent  a  near  free-molecular  flow  with  Kn  =4.81, 
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H  =  0.1  //m  and  AR=10.  Figure  1 1  shows  that  doubling  the  size  of  the  edge-length,  while  keeping  it 

smaller  than  ,  has  no  effect  on  centerline  flow  properties.  Case  4  and  Case  GS-4  correspond  to  a 

transitional  flow  with  Kn^  —  0.481,  //  =  0.5 //m  and  AR=1.  A  comparison  between  Case  4  and 

GS-4  flow  properties  in  Figure  1 1  shows  that  the  doubling  of  the  edge-length  has  no  effect  on  the  cen¬ 
terline  flow'  properties.  These  results  obtained  for  the  most  restrictive  cases,  provide  evidence  that  the 
choice  of  computational  parameters  used  in  the  nanoscale  simulations,  provide  sufficient  numerical 
resolution. 


1.5.3  Summary  of  Sim  ulation  Results 


1 . 5 .3 . 1  Effects  of  Nanochannel  Inlet  Size  and  Aspect  Ratio 


These  effects  were  examined  with  Case  1  to  Case  9  simulations  shown  in  Table  1.  The  inlet 
height  establishes  the  free  stream  Knudsen  number  and  characterizes  the  rarefaction  effects.  The  as¬ 
pect  ratio,  AR,  establishes  the  flow  development  due  to  particle-wall  interactions. 
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Figure  12.  Effects  of  the  aspeet  ratio  AR =///L  on  phase-space  (x,  y)%  (x%  Vy),  {vx,Vy)  obtained 

from  the  U3DSMC  simulation  of  a  supersonic  flow  into  an  //=0.1  pm  nanoehannel  (Case  1,  2  and  3).  The 
P  =0.1  atm  free  stream  is  to  the  left  of  the  1  pm  buffer  region  with  the  nanoehannel  inlet  indicated  by 
the  dotted  line. 
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The  phase  space  for  Cases  1 , 2  and  3  is  presented  in  Figure  12.  These  cases  simulate  the  smallest  na¬ 
nochannels  with  H  =  0.1 //m  and  correspond  to  the  near  free-molecular  flow  with  a  Kn  =4.81. 

The  (x,y)  phase  plot  shows  the  number  density  to  increase  inside  the  AR=1  nanochannel.  The 
AR=10,  and  AR=100  nanonchannels  show  an  increase  in  the  number  density  followed  by  a  decrease 
at  the  end  of  the  nanochannef  The  (.r,  vx )  phase  plot  shows  the  high-speed  free-stream  population  to 

permeate  the  entire  length  of  the  AR=1  nanochannel.  The  (x.v^)  plot  shows  that  inside  the  AR=1 

nanochannel  a  low-speed  population  develops  as  a  result  of  collisions  of  particles  with  the  walls.  In 
addition,  a  small  population  of  upstreaming  particles  exit  from  the  inlet  after  colliding  with  the  walls 

or  other  particles  in  the  density  enhancement  region.  The  {x.vj  plot  for  the  AR=10  and  AR=100 

cases  shows  that  the  high-speed  component  does  not  permeate  the  entire  length  of  the  nanochannel. 
Instead,  the  nanochannel  is  populated  with  near  zero-drift  particles,  due  to  diffuse  reflections  off  the 

walls.  The  (x,  v  )  phase  space  for  AR=10  and  AR=100  shows  that  the  spread  of  the  v  component  is 

reduced  towards  the  exit.  Both  phase  plots  exhibit  the  upstreaming  particles.  The  (x,v  )  phase  plots 

in  Figure  12  show  that  the  y-component  of  the  velocity  increases  inside  the  nanochannel  as  a  result  of 
the  collisions  in  the  density  enhancement  region.  As  the  AR  increases  from  1  tolOO,  the  velocity 

spread  decreases  with  increasing  length,  due  to  diffuse  reflections  off  the  walls.  The  (t;x,v  )  phase 

plots  in  Figure  12  corroborate  the  previous  results.  The  AR-1  nanochannel  shows  the  appearance  of  a 
symmetric,  free-streaming  population  and  a  secondary,  symmetric,  zero-drift  population  that  arise  due 
to  collisions  with  the  wall.  This  zero-drift  population  becomes  the  dominant  feature  of  the  flow  as  the 
AR  increases.  It  is  also  noticeable  that  the  spread,  that  is  the  temperature  of  the  low -drift  population, 
is  larger  than  the  free-stream  value  as  a  result  of  collisions  inside  the  nanochannel.  The  symmetry  of 
the  low-speed  population  is  indicative  of  gas-wall  equilibration  phenomena  inside  the  nanochannel. 

In  Figure  13  we  plot  the  sample-averaged,  macroscopic  fluid  properties  along  the  centerline  of 
the  nanochannel  and  the  buffer  region  They  include  the  number  density  from  eq.  (7),  the  average  ve¬ 
locity  V  from  eq.  (8),  and  the  translational  temperature  T  from  eq,  (11).  These  macroscopic  va¬ 
riables  are  sampled  in  Delaunay  cells,  shown  in  Figure  2,  that  occupy  a  considerable  fraction  of  the 
nanochannel  at  a  given  downstream  position.  With  a  few  sampling  cells  spanning  the  width  of  the  na¬ 
nochannel  it  is  not  feasible  to  obtain  a  meaningful  spatial  resolution  and  construct  contour  plots.  The 
number  density  shows  the  formation  of  the  density  enhancement  that  becomes  stronger  with  increas¬ 
ing  AR  for  all  Knudsen  numbers  considered.  The  density  decreases  inside  the  nanochannel  and  the 
trend  is  more  pronounced  with  the  larger  AR  cases. 
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Figure  13.  Centerline,  sample-averaged,  macroscopic  flow  properties  (FX,T,  n )  in  a  nanochannel  obtained 
from  U3DSMC  simulations  (Cases  1-9).  The  P ^  =  0.1  atm  free  stream  is  to  the  left  of  the  1  pm  buffer  re¬ 
gion  with  the  nanochannel  inlet  indicated  by  the  dotted  line. 


The  density  enhancement  is  explained  physically  with  the  phase-space  plots  discussed  earlier.  The 
average  velocity  decreases  near  the  inlet  of  the  nanochannel  for  all  cases  considered.  The  temperature 
shows  an  enhancement  in  the  inlet  region  and  it  is  associated  with  the  presence  of  the  high-speed  and 
low-speed  particle  populations  shown  in  the  phase  plots  of  Figure  12.  The  temperature  after  the  inlet 
region  shows  a  sharp  decrease  and  remains  almost  constant  for  the  AR=10  and  AR=100  nanochannels 

indicating  equilibrium  with  the  wall  temperature.  As  Table  1  shows  in  cells  with  nD  ^  the  num- 
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ber  of  computational  particles  is  —  10  —  40  and  in  cells  with  nD  >  the  value  of  (Nd}  can 

be  even  larger.  The  numerical  fractional  errors  in  the  sample-averaged  density,  velocity  and  tempera¬ 
ture  are  expected  to  be  below  the  range  given  in  Figure  5. 

For  all  cases  considered  the  density  maximum  is  located  near  the  inlet  but  inside  the  nanochan¬ 
nel.  The  temperature  maximum  is  located  before  or  at  the  inlet  of  the  nanochannel.  This  behavior 
where  density  and  temperature  enhancement  do  not  coincide  has  been  observed  in  rarefied  high-speed 
flows  into  macroscale  tubes  Gatsonis  et  al .  (1997)24  as  well  microchannels  (Oh  et  al 1997;25  Liou 
and  Fang,  200 1 ; ! <  Le  and  Hassan,  2007. r)  The  nanochannel  acts  a  solid  body  to  the  incoming  flow 
and  reflected  particles  from  the  wall  generate  this  perturbation  region  which  is  similar  to  a  diffuse 
shock  formed  in  front  of  a  solid  body  in  a  high-speed  rarefied  gas  flow  (Wilmoth  et  al.,  1 992)  29. 

In  Figure  14  the  centerline  pressure  from  eq.  (12)  and  Mach  number  are  plotted  for  Cases  1-9. 
The  AR=1  nanochannels  show  that  the  Mach  number  is  supersonic  at  the  inlet  and  remains  supersonic 
for  the  three  Knudsen  numbers  considered.  The  centerline  pressure  increases  from  its  free  stream 
value  and  decreases  slightly  towards  the  outlet.  The  nanochannels  with  AR=10  and  AR=100  result  in 
a  subsonic  inlet  Mach  number  and  a  sharp  increase  in  pressure  near  the  inlet.  The  Mach  number,  sub¬ 
sequently,  increases  and  reaches  the  sonic  point  at  the  outlet  for  all  cases  considered.  For  the  super¬ 
sonic  inlet  (Case  1,  4  and  7)  in  the  AR=1  nanochannels,  the  Mach  number  remains  supersonic  at  the 
outlet.  For  cases  where  the  inlet  Mach  number  becomes  subsonic  the  flow  becomes  sonic  at  the  exit. 
The  pressure  in  Figure  14  shows  a  sharp  increase  in  the  density  enhancement  region  near  the  inlet, 
followed  by  a  decrease  towards  the  end  of  the  channel.  Figure  13  and  Figure  14  show  the  1-pm  buf¬ 
fer  using  the  same  A-axis  size  for  simulation  domains  ranging  from  l.l  pm  (Case  1)  to  101  pm  (Case 
9).  Therefore,  for  some  of  the  AR-10  and  AR=100  cases  the  free  stream  conditions  cannot  be  de¬ 
picted  clearly.  While  there  are  variations  in  the  flow  parameters  in  the  buffer  region,  the  free  stream 
conditions  are  sustained  at  the  entry  of  the  l -pm  buffer  region  for  all  cases  considered. 

1 .5.3.2  The  Effects  of  Free  Stream  Speed  Ratio  and  Back  Pressure 

The  free  stream  speed  ratio  can  have  profound  effects  on  the  flow  characteristics  inside  a  nanochan¬ 
nel.  Simulations  for  H  =  0.5/i  and  AR=10  nanochannels  were  performed  with  free  stream  speed  ra¬ 
tios  of  S m  =  2  (Case  10)  ,  S ^  —  5  (Case  5)  and  —  10  (Case  1 1).  The  effects  of  backpressure 
were  examined  in  simulations  Case  4,  4a,  4b,  5,  5a,  5b  These  simulations  were  performed  with 
H  =  0.5//  and  AR=1  and  AR=10  nanochannels  with  free  stream  speed  ratio  of  S ^  =  5  and 

Px  =  10. 13  kPa.  Discussion  of  the  results  can  be  found  in  Al-Khouz  (2009)/ 
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Figure  14.  Centerline  sample-averaged  Mach  number  and  pressure  in  a  nanochannel  obtained  from 
U3DSMC  simulations  (Cases  1-9).  The  Px  =  101,320  pa,  =  5.97  free  stream  is  to  the  left  of  the  1 
pm  buffer  region  with  the  nanochanncl  inlet  indicated  by  the  dotted  line. 


1.5.4  Validation  of  U3DMC  at  the  Nanoscale 

The  heat  flux  to  the  wall  obtained  from  U3DSMC  through  eq.  (17)  is  plotted  in  Figure  15.  With 
the  exception  of  the  AR=1  nanochannels  the  heat  loss  is  largest  in  the  inlet  region.  A  comparison  be¬ 
tween  the  numerical  and  theoretical  values  of  heat  loss  can  be  obtained  for  Cases  1,2,3  that  corres¬ 
pond  to  near-free  molecular  flows  at  K n  =  4.51 .  The  free  molecular,  total  translational  and  inter¬ 
nal  energy  flux  (or  heat  transfer  rate)  from  a  free  stream  with  p<  T*S  to  a  flat  plate  with  surface 
temperature  T  aligned  with  the  flow  is. 
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where  £  is  the  fraction  of  molecules  that  is  reflected  specularly  (Bird,  2000)^.  Substituting  in  eq.  (26) 
values  for  p,T,S  from  the  U3DSCM  results  sampled  at  the  inlet  of  the  nanochannel  we  obtain  theo¬ 
retical  values  of  heat  flux.  For  Case  1,  =  —52.4  x  106  W/m“  with  the  numerical  q  =  —31.3  x  106 

W/irT.  For  Case  2,  q =  —93.7  x  106  W/rrf  and  q  =  —62.5  x  106  W/irf .  For  Case  3, 

q A  =  —42.5  x  106  and  q  =  -19.8  x  106  W/m  .  The  numerical  values  are,  therefore,  within  the  or¬ 
der  of  magnitude  of  the  theoretical  estimates  of  heat  loss. 
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Figure  15.  Heat  flux  to  the  wall  obtained  from  the  U3DSMC  simulation  of  a  supersonic  flow  into  a  nano¬ 
channel.  The  free  stream  is  at  A^O  pm  and  the  nanochannel  inlet  at  X=\  pm  (Case  1-9). 


The  outgoing  mass  flow  rate  at  the  outlet  from  U3DSMC  was  also  compared  with  estimates  ob¬ 
tained  using  the  semi-analytical  theory  developed  for  the  free-molecular  regime  by  Hughes  and  de 
Leeuw  (1965)21.  For  a  cylindrical  tube  with  aspect  ratio  D  =  d  /  L,  a  drifting  Maxwellian  with  a 
speed  ratio  at  zero  angle  of  attack  with  the  tube  axis,  will  have  a  flux  of  particles  exiting  at  the  outlet 
into  a  background  of  specified  pressure  given  by 


N(Sx,D.0) 


0 


K(<t>,  D)F(<t>,  ,  0)d(f> 


(27) 
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In  the  above,  Cm(x  is  the  most  probable  speed.  The  quantity  K(<p,D)  is  the  ratio  between  the  number 
flux  of  molecules  exiting  at  the  outlet  of  the  tube,  and  the  number  flux  of  particles  that  enter  the  tube 
due  to  a  beamlet  with  an  angle  0  with  the  tube  axis.  Hughes  and  de  Leeuw  evaluate  K(&  D)  using 
the  Clausing’s  probability  function  for  particle  transmission  after  wall  reflections.  The  flux  term  in  eq. 
(27)  that  strikes  a  surface  perpendicular  to  the  tube  axis  at  an  angle  between  0  and  d0 ,  is, 


F(0, 5,0)  =  sin  0  cos  0[(1  +  52  cos2  0)cxp(— 52)  + 

7 r1/25  cos0cxp(— 52  sin2  0)(|  -f  52  cos2  0)(1  +  erfS  cos0)] 


(28) 


The  expressions  in  eq.  (27)  have  been  integrated  numerically  and  evaluated  using  input  parameters 
from  the  simulations  for  Case  NO.  Comparisons  are  shown  in  Figure  16.  For  AR=1  and  AR=10,  the 
U3DSMC  predictions  are  very  close  to  the  theoretical.  For  the  long  AR=100  nanochannels,  the  theo¬ 
retical  predictions  are  larger  than  the  U3DSMC.  The  theory  does  not  account  for  the  external  flow  ef¬ 
fects  nor  the  internal  structure  of  the  flow  field  as  evident  in  the  U3DSCM  results  discussed  earlier. 
In  the  case  of  short  tubes,  wall  collisions  are  not  sufficient  to  modify  the  free-stream  and  the  theoreti¬ 
cal  results  provide  the  outgoing  flux. 
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Figure  16.  Outgoing  mass  flow  rate  at  the  outlet  of  a  nanochannel  for  various  aspect  ratios.  Comparison 
between  U3DSMC  simulations  (Case  N9)  and  numerical  results  obtained  from  the  Hughes  and  de  Teeuw 
,21(  theory. 


1 .6  U3DSMC  Simulation  of  Atmospheric-Pressure  Supersonic  Flows  Into  Nanochannels 

We  performed  simulations  of  supersonic  atmospheric  flows  into  nanochannels  using  the 
U3DSMC  method.  Figure  8  shows  the  simulation  domain  as  well  as  the  parameters  used  to  define  the 
nanochannel  and  the  buffer  region  geometries.  The  incoming  free  stream  of  nitrogen  N2  has 
n  =2.G9xl025in  \  T  =273K,  corresponding  to  P  =latm  and  A  =48.1  nm.  The  free 

stream  velocity  was  varied  in  the  simulations  so  that  the  speed  ratio  spans  the  range  of  the  super- 
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sonic  (S  ~  2)  to  hypersonic  (5^  =10)  flow  regimes.  The  variation  of  the  inlet  height  results  in 

Kn that  range  from  0.0962  to  0.481  covering  the  slip  to  transitional  regimes.  The  parameters  con¬ 
sidered  in  these  simulations  are  shown  in  Table  3.  The  simulations  considered  vacuum  back  pres¬ 
sure. 


Table  3.  Physical  and  computational  parameters  for  U3DSMC  simulations  of  atmospheric-pressure  ni¬ 
trogen  flows  into  a  nanochanncl. 


Ca 

se 

H 

(nm) 

LjH 

5 

oc 

M 

oc 

K 

(kPa) 

Kn 

oc 

Ls 

(nm) 

(nm) 

MU1 

Fs 

(*.) 
t  ■’(»„)] 

1 

100 

1 

5 

5.97 

0 

0.481 

1000 

1578 

8.62 

[0.793] 

12 

13.47 

[3.21] 

2 

100 

10 

5 

5.97 

0 

0.481 

1000 

2793 

8.45 

[0.775] 

12 

14.34 

[3.11] 

3 

500 

1 

5 

5.97 

0 

0.0962 

1000 

57407 

1 1.29 
[1.21] 

12 

14.21 

[3.07] 

4 

500 

10 

5 

5.97 

0 

0.0962 

1000 

57879 

18.08 

[2.15] 

30 

22.8 

[5.77] 

5 

100 

10 

2 

2.39 

1 1.9 

0 

0.0962 

1000 

2793 

8.45 

[0.775] 

8.45 

12 

14.34 

[3.11] 

14.34 

6 

100 

10 

10 

5 

0 

0.0962 

1000 

2793 

[0.775] 

12 

[3.11] 

The  phase  space  plots  for  Case  1  and  2  are  presented  in  Figure  17.  These  cases  simulate  the  narrowest 
nanochannels  with  H  =  0.1  /mi  and  correspond  to  the  transitional  flow  with  a  Kn ^  =0.481.  The 

(x,y)  phase  plot  shows  that  there  is  an  enhancement  in  the  number  density  inside  the  AR=1  nano¬ 
channel.  The  long  AR=10  nanonchannel  shows  an  increase  in  the  number  density  followed  by  a  de¬ 
crease  at  the  end  of  the  nanochannel.  The  (x,vx)  phase  plot  shows  the  high-speed  free-stream  popula¬ 
tion  to  permeate  the  entire  length  of  the  AR=1  nanochannel.  The  (x,  v  )  plot  shows  a  low-speed 

population  of  particles  that  develops  inside  the  AR=1  nanochannel  as  a  result  of  collisions  of  particles 
with  the  walls.  In  addition,  a  small  population  of  up  streaming  particles  exit  the  inlet  after  colliding 
with  the  walls  or  other  particles  in  the  density  enhancement  region.  The  (x,vx)  plot  for  the  AR=10 

case  show  that  the  high-speed  component  does  not  permeate  the  entire  length  of  the  nanochannel.  In¬ 
stead,  the  nanochannel  is  populated  with  a  near  zero-drift  particles  due  to  diffuse  reflections  of  the 

walls.  The  (x,vx)  phase  space  shows  that  the  spread  of  the  vx  component  is  reduced  towards  the  exit. 
The  (x,v  )  phase  plots  show  that  the  y-component  of  the  velocity  increases  inside  the  nanochannel  as 
a  result  of  the  collisions  in  the  density  enhancement  region.  As  the  AR  increases  the  velocity  spread 
decreases  with  increasing  length,  due  to  diffuse  reflections  of  the  walls.  The  (vx,vy)  phase  plots  cor- 
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roborates  the  analysis  of  the  space-velocity  phase  plots.  The  AR=1  nanochannel  show  the  appearance 
of  a  symmetric  free-streaming  population  and  a  secondary,  symmetric,  zero-drift  population  that  arise 
due  to  collisions  with  the  wall.  This  zero-drift  population  becomes  the  dominant  feature  of  the  flow 


as  the  AR  increases.  It  is  also  noticeable  that  the  spread  in  v  (which  is  an  indication  of  the  tempera¬ 
ture)  for  the  low-drift  population  is  larger  than  the  free-stream  value  as  a  result  of  collisions  inside  the 
nanochannel.  The  symmetry  of  the  low-speed  population  is  indicative  of  gas-wall  equilibration  phe¬ 
nomena  inside  the  nanochannel.  The  complete  set  of  results  and  discussion  is  presented  in  Al-Khouz 


(2009). 
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Figure  17.  FTfects  of  AR  for //=0.1  pm  nanochannels.  Phase  plots  (x,y)*  (x,  v* ) ,  (x,  v )  and  (va,t;  )  for 

Case  1,  2. 


1  7Conclusions 

The  U3DSMC  simulations  in  this  work,  featured  a  small  number  of  real  particles  in  the  computa¬ 
tional  domain,  small  particle  weights,  and  few  computational  cells  across  the  width  of  the  nanochan- 
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nel.  The  U3DSMC  results  showed  that  the  supersonic  nanoflows  considered  exhibit  characteristics 
that  are  similar  to  previous  micron-  and  macroscale  flows  covering  similar  Knudsen  number  regimes, 
and  have  physical  characteristics  that  are  found  in  continuous  compressible,  viscous  channel  flows 
with  friction  and  heat  loss.  The  microscopic  analysis  adopted  in  our  work  used  the  phase-space  dis¬ 
tributions  and  provided  the  kinetic-based  explanation  of  the  obtained  macroscopic  flow  variables. 
Overall,  our  work  extended  the  applicability  of  the  U3DSMC  method  to  rarefied  supersonic  nanof¬ 
lows. 

2  SMOOTHED  DISSIPATIVE  PARTICLE  DYNAMICS 

We  review  below  the  second  major  approach  undertaken  during  this  research  involving  particle  simu¬ 
lation  of  fluid  flows  at  the  mesoscale.  The  Smooth  Dissipative  Particle  Dynamics  (SDPD)  method 
developed  by  Espanol  and  Revenga,  (2003)  '"  is  a  Lagrangian,  mesh-free,  particle  model  that  can  be 
applied  from  microscale  to  nanoscale  flow  domain,  as  shown  in  Figure  18.  SDPD  can  be  derived  as  a 
top-down  approach  from  SPH  with  the  introduction  of  a  fluctuation  term.  As  a  bottom-up  approach, 
SDPD  it  is  a  thermodynamically  consistent  version  of  the  original  DPD,  due  to  the  introduction  of  the 
energy  equation  in  the  model. 
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Figure  18.  Length  scales  and  computational  approaches. 


In  SDPD  the  fluid  particle  with  position  r  and  velocity  v  has  a  volume  defined  as 

1  JN 

—  =  (l  =  V'  W(r  ,  A)  with  r.  =1  r  —  r.  I 

y  t  v  17  ■  '  i]  1  «  ]  1 


(29) 


The  supporting  function  is  similar  to  the  one  found  in  SPH  One  choice  based  on  the  interaction 
length  A  is 
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The  density  of  the  i-th  SDPD  particle  is 

N 

Pi  =  mA  =  miYlw(ri)'h') 

i- 1 


(30) 
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The  mass 


Mi 

The  SDPD  particle  position  is 
dr  =  v  dt 

t  t 

The  SDPD  momentum  is  given  by 


(32) 


(33) 


(34) 


The  first  term  in  the  momentum  equation  provides  the  conservative  force.  The  pressure  of  an  SDPD 
particle  is  derived  from  the  equation  of  state 
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The  second  group  of  terms  in  momentum,  provides  the  dissipative  force  with  a  rotational  and  central 
component.  The  third  term  provides  the  random  force 


mdvt  = 


A  dW,'  +  n  -«r[dW 


(36) 


The  dW  is  a  matrix  of  independent  increments  of  the  Wiener  process.  Dissipative  and  random 

forces  are  related  in  SDPT  through  the  temperature  terms. 

The  entropy  equation  for  the  i-the  SDPD  particle  is 
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Details  for  the  coefficients  in  the  SDPD  model  equations  are  given  in  Espanol  and  Revenga  (2003), 
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2.1  SDPD  Computational  Implementation 

The  SDPD  model  has  been  implemented  in  a  serial  version  of  a  code  shown  in  Figure  19.  The 
implementation  has  the  capability  of  handling  complex  geometries.  The  integration  scheme  for  the 
equation  of  SDPD  particle  motion  (33)  is  the  modified  velocity-Verlet  scheme  (Groot  and  Warren, 
1997)  2  where  the  dissipative  force  component  is  calculated  in  two  steps  for  each  integration  cycle. 
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Periodic  and  solid  wall  boundary  conditions  have  been  implemented  into  the  code.  The  entropy  equa¬ 
tion  has  also  been  implemented  and  undergone  preliminary  validations. 


MAIN 


Module  1 


2.2  Solid  Wall  Boundary  Conditions 

The  implementation  of  solid  wall  boundary  conditions  in  DPD  methods  is  an  active  research  area 
(Fedosov  et  al,  2009).  We  developed  under  this  grant  an  SDPD  approach  to  represent  the  geometry 
of  a  solid  wall.  In  our  model  the  solid-wall  is  represented  by  SDPD  particles  with  positions.  The  sol¬ 
id-wall  SDPD  particle  velocity  is  chosen  to  be  opposite  in  sign  and  having  the  same  magnitude  as  the 
interacting  fluid  particle  velocity.  In  the  current  implementation  of  SDPD  the  solid-wall  particle  does 
not  affect  the  properties  of  the  neighboring  fluid  particle  but  does  affect  the  acting  force.  In  our  mod¬ 
el,  the  force  of  the  solid-wall  and  fluid  particles  is  “strong  enough”  to  avoid  penetration  of  fluid  par- 
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tides  inside  the  wall  without  the  need  of  arbitrary  constants  as  has  been  the  case  of  DPD  implementa¬ 
tions.  The  algorithmic  steps  are  summarized  below  and  the  geometrical  characteristics  are  shown  in 
Figure  20  . 

A  conservative  force  component  is  computed  first: 

S-l.  Compute  the  conservative  force  Fr  acting  on  particle  1  from  its  neighbor  fluid  par¬ 
ticles  contained  in  its  finite  support  h 
S-2.  Compute  the  1  particle-wall  unit  vector  e 

S-3.  Compute  the  projection  of  the  conservative  force  from  step  1)  on  direction  e 

S-4  If  the  projected  force  from  Fc  is  directed  towards  the  wall  then  apply  this  contribu¬ 

tion  to  the  i-th  particle,  F^1  =  ( e  •  F^  )  e 


A  dynamic  force  component  is  computed  next: 

S-l .  Compute  the  fluid  velocity  component  V*  along  direction  e 

S-2.  Assume  that  at  the  wall  V*  =  0 

S-3.  Compute  the  dynamic  component  for  the  total  force  as  function  of  the  module  of  the 


normal  velocity  component,  F^2 


777  A  VN 

— - 1 — e 

At 


S-4.  If  the  normal  velocity  component  is  directed  towards  the  wall  then  apply  this  contri¬ 
bution  to  the  i-th  particle 

S-5.  Force  on  the  i-th  particle  from  the  wall  is  Fc  —  FC1  +  FC2 

1  t  w  w 


Figure  20.  Solid  wall  implementation  in  SDPD. 


Compute  a  dissipative  force  component  for  the  i-th  particle: 

S-l .  For  a  defined  distance  Calculate  within  a  defined  distance  from  the  wall  compute  the 
velocity  component  V7  normal  to  e  (tangential  to  the  wall) 

S-2.  Assume  that  at  the  wall  Vr  =  0 


S-3. 


Compute  the  force  due  to  V7  given  by  F^  =  —m 


AVT 

_ t_ 

’  At 


and  apply  to  the  particle. 
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2.3  Validation  of  the  SDPD  Method 


The  objective  of  this  simulation  was  to  validate  the  3D  SDPD  code  as  well  as  the  solid-wall 
boundary  condition  approach.  An  incompressible  flow  was  considered  in  a  microchannel  with  circu¬ 
lar  constant  section  of  radius  9  microns  and  length  of  20  microns.  The  fluid  simulated  is  air  with  vis¬ 
cosity  of  7/  —  10  6  Kg  /  ms ,  temperature  T  =  288  K  ,  density  of  0.7 6Kg  /  m 3  and  a  Reynolds  num¬ 
ber  of  50.  The  total  number  of  simulate  SDPD  fluid  particles  was  6561  which  implies  for  the  given 
fluid  a  coarsening  factor  Nm  ~  107  and  a  total  Knudsen  number  of  Kn  ~  10  z .  The  flow  was  driven 

by  imposing  an  external  force  acting  on  each  SDPD  fluid  particle.  Frozen  SDPD  particles,  with  as¬ 
signed  position,  were  used  to  represent  the  solid  wall  of  the  microchannel.  Periodic  conditions  were 
assigned  at  the  ends  of  the  channel. 

Figure  21  shows  the  axial  velocity  from  SDPD  along  with  the  analytical  solution  for  Poiseuille 
flow.  Results  show  very  good  agreement  between  the  two  profiles. 


Figure  21.  Axial  velocity  from  SDPD  and  comparison  with  analytical  Poiseuille  profile. 

A  common  issue  integrating  the  particle  equation  of  motion  for  mesoscopic  dynamics  is  the  ex¬ 
hibit  of  nonphysical  behavior  such  as  the  systematic  drift  in  temperature  (for  isothermal  cases).  In  or¬ 
der  to  control  the  quality  of  the  integration  scheme  the  averaged  kinetic  energy  was  monitored  as 
shown  in  Figure  22.  The  results  display  the  good  convergence  and  small  level  of  fluctuations  of  our 
SDPD  implementation.  The  density  across  the  microchannel  shown  in  Figure  22  has  a  fractional  fluc¬ 
tuation  of  less  than  0.22%  and  demonstrates  the  ability  of  the  solid-wall  boundary  condition  devel¬ 
oped  in  this  work. 
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Figure  22.  Average  kinetic  energy  as  a  function  of  time  and  density  across  the  channel  from  the 

SDPD  simulation. 


A  second  simulation  example  involved  a  fluid  in  a  periodic  domain  with  given  temperatures  and 
random  velocities.  The  entropy  equation  (37)was  included  in  the  simulation  with  the  heat  transfer 
term  shown  in  eq.  (37).  The  simulation  was  run  for  about  500  timesteps. 

Figure  22  shows  the  entropy  and  average  SDPD  particle  temperature  in  the  domain.  The  in¬ 
crease  in  entropy  and  the  average  temperature  is  due  to  the  irreversibility  of  the  heat  transfer  term  in 
the  entropy  equation.  The  simulation  shows  that  equilibration  is  maintained  without  any  unphysical 
increase  in  temperature. 
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Figure  23.  Entropy  and  average  temperature  from  SDPD  simulation. 


2.4  Conclusions 

We  implemented  successfully  the  SDPD  method,  a  thermodynamically  consistent  DPD  model.  We 
developed  a  boundary  condition  methodology  that  does  not  require  arbitrary  constants.  The  method 
and  code  have  been  validated  successfully  for  isothermal  pressure  driven  flows.  Simulations  of  non- 
isothermal  flows  exhibited  good  convergence  properties. 
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